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The  APEX  method  is  a non-iterative,  single  frame,  direct  blind  deconvolution 
technique  that  can  sharpen  certain  kinds  of  high  resolution  images  in  quasi 
real-time.  The  method  is  predicated  on  a restricted  class  of  blurs,  in  the 
form  of  2D  radially  symmetric,  bell-shaped,  heavy-tailed  probability  density 
functions.  Not  all  images  can  be  usefully  enhanced  with  the  APEX  method. 
However,  the  method  is  found  effective  on  a broad  class  of  galaxy  images, 
including  color  Hubble  Space  Telescope  Advanced  Camera  for  Surveys  (ACS) 
imagery.  APEX-detected  optical  transfer  functions  that  successfully  sharpen 
these  images  are  very  far  from  Gaussian,  and  of  a type  seldom  found  in  the 
imaging  literature.  Significantly  sharper  arid  visually  striking  reconstructions 
are  obtained,  with  sharpening  confirmed  by  the  tripling  or  quadrupling  of 
image  gradient  norms.  © 2005  Optical  Society  of  America 

OCIS  codes:  100.1830,  100.3020,  110.4850,  110.6770. 


1.  Introduction 

The  APEX  method  is  a non-iterative,  single  frame,  direct  blind  deconvolution  technique 
that  can  sharpen  certain  kinds  of  high  resolution  images  in  quasi  real-time.  The  method 
operates  in  Fourier  transform  space  via  FFT  algorithms.  Typically,  1024  x 1024  pixel  images 
can  be  processed  in  seconds  on  current  desktop  computers.  The  method  has  been  applied 
successfully  in  diverse  imaging  contexts,  including  airborne  reconnaissance,  MRI  and  PET 
brain  scans,  and  scanning  electron  microscopy.1-3  However,  not  all  images  can  be  usefully 
enhanced  with  the  APEX  method.  The  present  paper  explores  the  possible  application  of 
this  technique  to  astronomical  data,  including  Hubble  Space  Telescope  (HST)  imagery.  In 
Figure  1,  a familiar  earthbound  setting  illustrates  the  type  of  improvement  that  is  sometimes 
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Fig.  1.  APEX  blind  deconvolution  of  English  village  image,  (a)  Original  512  x 
512  8-bit  image,  fb)  APEX  processed  image.  Zooming  on  selected  parts  of 
sharpened  image  (b)  reveals  buildings  in  the  distance,2  and  other  significant 
information  not  easily  detectable  in  image  (a). 


Fig.  2.  Hubble  Space  Telescope  image  of  M100  galaxy,  before  and  after  imple- 
mentation of  1993  corrective  optics  package.  APEX  method  applied  to  image 

o 

(a)  would  fail  to  detect  flawed  optics  psf,  and  would  be  unable  to  produce 
useful  approximation  to  image  (b). 


possible  with  the  APEX  method.  In  that  example,  zooming  on  selected  parts  of  the  APEX- 
enhanced  image  1(b)  reveals  buildings  in  the  distance,2  Holstein  cows  grazing  in  the  meadow, 
and  numerous  other  fine-scale  details  not  readily  apparent  in  the  original  image  1(a). 

In  recent  years,  much  excellent  work  has  been  done  in  the  area  of  blind  deconvolution 
of  astronomical  data.4-11  Many  of  these  methods  aim  primarily  at  undoing  the  distorting 
effects  of  atmospheric  turbulence  in  short-exposure,  ground-based  observations.  Multiframe 
algorithms,  typically  involving  several  hundred  short-exposure  images  of  the  same  object, 
appear  to  be  particularly  effective.  An  interesting  example  of  multiframe  blind  deconvolution, 
applied  to  ground-based  surveillance  of  space  objects,  is  given  in  Ref.  9. 

APEX  processing  is  generally  not  useful  in  such  short-exposure  applications,  and  the 
method  would  probably  be  incapable  of  reproducing  the  results  in  Ref.  9.  In  a similar  vein, 
consider  the  severely  blurred  early  Hubble  Space  Telescope  imagery  caused  by  manufacturing 
flaws  in  the  primary  mirror.  The  much  improved  imagery  following  the  1993  implementation 
of  corrective  optics  is  best  illustrated  with  the  M100  galaxy  images  in  Figure  2.  Here,  if 
APEX  processing  were  to  be  applied  to  the  blurred  image  in  Figure  2(a),  the  method  would 
fail  to  identify  the  flawed  optics  point  spread  function  from  the  data  in  2(a),  and  it  would 
be  unable  to  produce  a useful  approximation  to  the  sharp  image  in  Figure  2(b). 

The  APEX  method  is  predicated  on  an  important  but  circumscribed  class  of  radially 
symmetric  shift- invariant  blurs,  one  that  generalizes  Gaussian  and  Lorentzian  distributions. 
This  is  the  class  G defined  in  Eq.  (4)  below.  That  class  does  not  include  the  more  complex 
point  spread  functions  that  characterize  the  examples  mentioned  in  the  preceding  paragraph. 
Rather,  the  APEX  method  aims  primarily  at  reconstructing  fine  scale  information  that  may 
have  been  smoothed  out  by  the  combined  effects  of  radially  symmetric  lens  aberrations,  long- 
exposure  turbulence  if  present,  and  additional  radially  symmetric  blurring,  originating  from 
diverse  electron  optical  devices  used  in  the  acquisition  and  recording  of  the  final  digitized 
image.  Presumably,  the  APEX  method  will  be  successful  on  a given  image,  only  to  the 
extent  that  a significant  portion  of  the  unknown  image  blur  can  be  well- approximated  by 
some  member  of  the  class  G. 

It  develops  that  APEX  processing  is  surprisingly  effective  on  a broad  class  of  galaxy 
images,  including  color  Hubble  Space  Telescope  imagery.  APEX-detected  optical  transfer 
functions  that  successfully  sharpen  these  images  are  very  far  from  Gaussian,  and  of  a type 
not  commonly  found  in  the  astronomical  imaging  literature.  Visually  striking  enhancements 
are  exhibited  in  sections  9 through  11.  The  degree  of  sharpening  in  these  images  can  be  quan- 
titatively assessed  by  comparing  image  gradient  norms  before  and  after  APEX  processing. 
Tripling  or  quadrupling  of  gradient  norms  is  commonly  realized.  In  the  case  of  Hubble  Tele- 
scope imagery,  the  APEX  method  can  enhance  Advanced  Camera  for  Surveys  (ACS)  images, 
in  addition  to  Wide  Field  and  Planetary  Camera  2 (WFPC2)  images.  An  interesting  new 
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method  of  assessing  image  sharpness,  based  on  measuring  image  Lipschitz  exponents ,12  can 
also  be  fruitfully  applied  to  the  present  class  of  images.  However,  due  to  space  limitations, 
a full  discussion  of  that  technique  must  be  deferred  to  a future  report. 


2.  Heavy-tailed  Levy  point  spread  functions 


Important  empirical  work  has  identified  the  general  functional  form  of  the  optical  transfer 
functions  (otf)  in  a very  wide  variety  of  electron  optical  imaging  devices,13,14  including 
phosphor  screens  and  some  types  of  photographic  film.  Define  the  2D  Fourier  transform 
tF{h}  of  a function  h(x,y)  by 

tF{h}  = h(£,  r?)  = / h(x,  y)exp{— 2iri(^x  + r]y)}dxdy.  (1) 

Jr2 

When  h{x,y)  is  a point  spread  function  (psf),  it  is  non-negative  and  integrates  to  unity. 
Such  a function  corresponds  to  a probability  density  function.  The  optical  transfer  function 
h(£,rj)  corresponds  to  the  characteristic  function  of  that  density.  Apparently,  most  electron 
optic  imaging  devices  have  otfs  that  can  be  expressed  by 

h(^rj)  = exp{—a(£2  + rj2)3}.  a > 0,  0 < (3  < 1.  (2) 


where  the  constants  q and  3 depend  on  the  particular  device.13- 14  The  corresponding  densities 
h(x , y)  are  bell-shaped  surfaces  in  physical  x , y space,  and  belong  to  the  class  of  radially 
symmetric  Levy  stable  laws. 15,16  The  constant  a > 0 in  Eq.  (2)  controls  the  width  of  the 
density  h(x,y ),  and  h(x,y)  approaches  the  Dirac  5-function  as  a — ► 0.  The  constant  (3  is 
called  the  Levy  exponent.  The  case  (3  = 1 in  Eq.  (2)  corresponds  to  the  Gaussian  distribution, 
while  the  case  (3  — 1/2  corresponds  to  the  2D  Lorentzian  density 

a 


h(x,y)  = - 


2x(x2  + y2  + a2)3'2' 


(x,  y ) e R‘ 


(3) 


For  other  values  of  3 in  Eq.  (2),  the  corresponding  density  h(x,  y ) is  not  known  in  closed  form 
in  the  physical  variables  x,  y.  In  the  Gaussian  case  (3  = 1,  h(x,  y)  has  exponentially  decaying 
slim  tails  and  finite  variance.  However,  for  0 < [3  < 1,  h(x.  y)  has  infinite  variance,10-19  with 
heavy  tails  that  decay  like  a power  of  1/r,  where  r = (x2  + y2)12  ■ 

The  expression  in  Eq.  (2)  can  be  used  to  describe  other  important  types  of  blurs.  The  otf 
for  long-exposure  turbulence  blurring  is  given  by  Eq.  (2)  with  3 = 5/6  and  a determined 
by  atmospheric  conditions.20  The  analytically  known  diffraction-limited  otf  for  a perfect  lens 
can  be  approximated  over  a wide  frequency  range  by  Eq.  (2),  with  3 = 3/4  and  a a properly 
chosen  function  of  the  cutoff  frequency.21  Otf  data  for  56  different  kinds  of  photographic  film 
have  also  been  analyzed.22  Good  agreement  is  found  when  these  data  are  fitted  with  Eq.  (2). 
and  the  pairs  (a,  (3)  characterizing  each  of  these  56  otfs  are  identified.  It  is  found  that  36 
types  of  film  have  otfs  where  1/2  < (3  < 1.  The  remaining  20  types  have  values  of  3 in  the 
range  0.265  < (3  < 0.475.  Corresponding  psfs  are  very  far  from  Gaussian. 


3.  Generalized  Central  Limit  Theorem  and  the  APEX  method 

The  classical  Central  Limit  Theorem  considers  the  limiting  probability  distribution  of  nor- 
malized sums  of  large  numbers  of  independent  random  variables  with  finite  variance , and  it 
asserts  that  that  limit  is  always  a Gaussian  distribution.1  ’ In  fact,  Gaussians  are  often  used 
to  fit  empirically  obtained  bell-shaped  data,  and  this  choice  is  usually  justified  on  the  basis 
of  that  theorem.  For  an  example  of  just  such  an  approach  applied  to  electron  optics  point 
spread  functions,  see  Ref.  23. 

In  recent  years,  with  the  advent  of  more  sophisticated  measurement  methods,  numerous 
physical  situations  have  been  uncovered  where  Gaussians  provide  inadequate  descriptions 
of  observed  bell-shaped  data,  because  legitimate  heavy-tailed  behavior  cannot  be  accomo- 
dated.1' l!l  21  A particularly  instructive  discussion  of  the  need  to  consider  non-Gaussian  dis- 
tributions is  contained  in  a recent  article  on  high  energy  particle  physics.25  It  is  now  gener- 
ally recognized  that  such  heavy- tailed  data  reflect  underlying  random  processes  with  infinite 
variance . and  that  such  processes  are  pervasive  in  nature.26  Evidently,  the  empirical  work 
reported  in  Refs.  13,  14.  and  22,  is  merely  one  instance  of  a recurring  pattern. 

The  Generalized  Central  Limit  Theorem  considers  normalized  sums  of  independent,  iden- 
tically distributed  random  variables,  with  variances  that  need  not  be  finite.  According  to 
that  theorem,  the  limit  of  any  such  sum,  if  it  exists,  must  be  a Levy  stable  law.15,18  Note 
that  while  the  class  of  stable  laws  includes  more  complex  asymmetric  specimens,  this  paper 
restricts  attention  to  the  radially  symmetric  case  through  Eq.  (2) 

In  some  applications,  several  electron  optical  devices  may  be  cascaded  together  and  used 
to  image  objects  through  a distorting  medium  such  as  the  atmosphere.  The  overall  lumped 
psf  is  then  the  convolution  product  of  the  individual  component  psfs,  so  that 

j 

h^.r/)  = exp{-Y,®i(£2  + ?/2)^},  W > 0,  0 < ft  < 1.  (4) 

i= 1 

The  general  functional  form  given  in  Eq.  (4)  may  also  be  used  to  best-fit  a large  class 
of  empirically  determined  optical  transfer  functions,  by  suitable  choices  of  the  parameters 
cq,  fii,  and  J. 

We  define  the  class  G of  blurring  kernels  to  be  the  class  of  all  psfs  h(x , y)  whose  Fourier 
transforms  satisfy  Eq.  (4).  We  shall  be  interested  in  image  deblurring  problems 

Hf  = h(x  — u,y  — v)f(u , v)dudv  = h(x,  y ) 0 f{x,  y)  = g(x,  y),  (5) 

■JR2 

where  g(x.y)  is  the  recorded  blurred  image,  f(x.y)  is  the  desired  unblurred  image,  and 
h(x.y)  is  a known  point  spread  function  in  class  G.  The  blurred  image  g(x,y)  includes 
noise,  which  is  viewed  as  a separate  additional  degradation, 

g(x,  y ) - ge(. r,  y)  + n(x,  y).  (6) 
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Here,  ge(x,y)  is  the  blurred  image  that  would  have  been  recorded  in  the  absence  of  noise, 
and  n(x,  y)  represents  the  cumulative  effects  of  all  errors  affecting  final  acquisition  of  the 
digitized  array  g(x,y).  The  unique  solution  of  Eq.  (5)  when  the  right  hand  side  is  ge(x,y ), 
is  the  exact  sharp  image  denoted  by  fe(x,y).  Thus 

h(x,  y)  <g>  fe{x,  y)  = ge(x,  y).  (7) 

Class  G otfs  are  non-negative.  This  is  not  generally  the  case  with  characteristic  functions. 
With  class  G otfs  we  may  define  fractional  powers  H\  0 < t < 1.  of  the  convolution  integral 
operator  H in  Eq.  (5)  as  follows 

= 0 < t < 1.  (8) 

Class  G psf's  are  intimately  related  to  diffusion  processes,  in  that  u(x,  y , t)  = f is  the 
solution  at  time  t of  a generalized  diffusion  equation,  (see  Eq.  (13)  below). 

These  considerations  underlie  the  APEX  blind  deconvolution  approach,  which  stipulates  at 
the  outset  that  the  blurring  is  isoplanatic,  and  that  the  lumped  total  system  optical  transfer 
function  can  be  well  approximated  by  Eq.  (4).  The  APEX  method  is  based  on  detecting 
such  Levy  stable  psfs  by  appropriate  Fourier  analysis  of  the  blurred  image  data.  As  discussed 
more  fully  below,  detected  representative  values  for  the  constants  cq  and  ft  in  Eq.  (4)  are 
used  to  construct  a candidate  otf.  This  is  then  used  in  the  SECB  deconvolution  method, 
implemented  as  a time-reversed  diffusion  equation.  By  marching  backwards  in  time,  one  can 
visually  monitor  the  deconvolution  process  as  it  unfolds,  examine  accompanying  diagnostic 
information,  and,  if  necessary,  choose  to  terminate  that  process  prior  to  completion.  Early 
termination  is  equivalent  to  interactive  readjustment  of  the  initial  candidate  otf. 

4.  Images  and  their  Fourier  transforms 

The  Fourier  transform  is  the  primary  computational  tool  used  in  this  paper.  We  deal  exclu- 
sively with  square  images  g(x,  y)  of  size  2 N x 2N  pixels.  In  order  to  render  mathematical 
formulae  more  transparent,  we  use  the  same  notation.  <?(ft  ?/),  for  both  discrete  and  continu- 
ous Fourier  transforms.  In  the  discrete  FFT  case,  the  frequencies  £ and  ij  are  understood  to 
be  integer-valued  and  to  range  from  — N to  N.  Likewise,  g(x,  y)  denotes  both  discrete  and 
continuous  images.  In  the  discrete  case,  the  variables  x , y are  measured  in  pixels  and  range 
from  1 to  2N . 

Given  the  Levy  pairs  (a*,  ft),  i = 1,  J,  where  cq  > 0,  0 < ft  < 1,  the  corresponding 
discrete  class  G otf  is  the  2Ar  x 2AT  array  h(£,rj)  where,  with  integer  ft  i) 

He  V)  = exp{-J2a i(?2  + iff1}- 

i=l 


N <£,T)<  Ar. 


(9) 


Fig.  3.  Fourier  behavior  in  1024  x 1024  image  of  spiral  galaxy  M101  is  typical 
of  a large  class  of  astronomical  images.  Above  image  was  taken  by  Jacoby,  Bo- 
hannan  and  Hanna,  Kitt  Peak  National  Observatory,  (NOAO/AURA/NSF). 
(a)  In  |p*(f.0)|  on  |£|  < 500  for  M101  image.  While  local  behavior  is  highly 
oscillatory,  global  behavior  is  generally  monotone  decreasing  and  convex,  (b) 
Least  squares  fit  of  In  |<?*(£,  0)|  with  u(£)  = — q|£|2J  — A , with  A = 3.85.  Fit 
develops  cusp  at  £ = 0 and  returns  a = 0.385,  = 0.206.  As  explained  in 

section  8,  such  trial  least  squares  fits,  using  different  values  of  A > 0,  are  basic 
to  the  APEX  method. 
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NATURAL  LOG  OF  ABS  VALUE  OF  FFT 


In  this  paper,  typical  parameter  values  might  be  TV  = 512,  J — 1,  a = 0.2,  j3  = 0.2.  Such 
otf  arrays  are  used  to  construct  the  SECB  deblurred  image,  as  in  Eq.  (14)  below. 

Given  an  image  g(x,  y),  the  natural  logarithm  of  the  absolute  value  of  its  Fourier  transform, 
In  |g(G7?)l>  wiH  play  a crucial  role.  This  logarithm  is  well-defined  except  where  #(£,77)  — 0. 
At  any  such  zero,  we  simply  redefine  g to  be  the  machine  epsilon.  In  practice,  exact  zeroes 
of  y(£ , 77)  are  seldom  encountered  due  to  system  noise. 

The  qualitative  behavior  in  Fourier  space  of  a large  class  of  astronomical  images  is  of 
interest.  Let  fe(x,y)  be  an  exact  sharp  image  as  in  (7).  Since  fe(x,y)  > 0 

\fe(^g)\<  I fetx,y)dxdy  = fe{ 0,0)  =7  >0.  (10) 

JR2 

Also,  since  ge(x,y)  = h(xyy)  0 fe(x,y)  and  h(x,y)  is  a probability  density, 

<7e(0, 0)  = / ge(x,y)dxdy=  fe(x,y)dxdy 

Jr 2 Jr2 

= fe( 0,0)  = 7 >0.  (11) 

Using  7 as  a normalizing  constant,  we  may  normalize  any  Fourier  transform  quantity  g(£,  //) 
by  dividing  by  7.  Let 

<?*(£/'?)  = </(£■>?)/ 7,  (12) 

denote  the  normalized  quantity.  The  function  |/e  ($,7/)|  is  highly  oscillatory,  and  0 < 
|/e  | A L Since  fe(x,y ) is  real,  its  Fourier  transform  is  conjugate  symmetric.  Therefore,  the 
function  \fe  ($,  //)  | is  symmetric  about  the  origin  on  any  line  through  the  origin  in  the  (£,77) 
plane.  The  same  is  true  for  the  blurred  image  data  |p*(£,77)|. 

For  any  2 TV  x 2 TV  image  g(x,  y ),  the  discrete  FFT  y(£,  77)  is  a 2 TV  x 2 N array  of  complex 
numbers.  The  frequencies  are  integers  lying  between  —N  and  TV , and  the  zero  frequency 
is  at  the  center  of  the  transform  array.  This  ordering  is  achieved  by  pre-multiplying  g(x,  y) 
by  ( — l)x+y.  We  shall  be  interested  in  the  values  of  such  transforms  along  single  lines  through 
the  origin  in  the  discrete  (£,77)  plane.  The  discrete  transforms  |<?*(f,0)|,  and  \g*(0,rj)\  are 
immediately  available.  Image  rotation  may  be  used  to  obtain  transforms  along  other  direc- 
tions. All  1-D  Fourier  plots  shown  in  this  paper  are  taken  along  the  axis  77  = 0 in  the  (£,  77) 
plane,  as  is  the  case  in  Figure  3.  In  these  plots,  the  zero  frequency  is  at  the  center  of  the 
horizontal  axis,  and  the  graphs  are  necessarily  symmetric  about  the  vertical  line  £ = 0. 

The  class  of  astronomical  images  g(x,y)  considered  in  the  present  paper  can  be  described  in 
terms  of  the  behavior  of  In  \g*(£,  77)  | along  single  lines  through  the  origin  in  the  (£,  77)  plane. 
While  local  behavior  is  highly  oscillatory,  global  behavior  is  generally  monotone  decreasing 
and  convex  on  £ > 0.  This  is  shown  in  Figure  3(a)  for  a typical  galaxy  image  along  the  line 
77  = 0,  and  similar  behavior  is  found  along  other  lines  through  the  origin  in  the  (£,77)  plane. 
A least,  squares  fit  to  the  oscillatory  trace  in  Figure  3(a)  with  a smooth  curve,  provides  a 
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good  representation  of  this  global  monotone  convexity  property  on  £ > 0.  (A  convex  function 
is  such  that  given  any  two  distinct  points  A and  B on  its  graph,  the  straight  line  segment 
joining  A and  B lies  above  the  graph.)  Many,  but  not  all,  astronomical  images  exhibit  similar 
globally  monotone  convex  Fourier  behavior.  Figure  3(b)  illustrates  the  type  of  least  squares 
fit  that  is  basic  to  the  APEX  method,  and  that  will  be  described  more  fully  in  section  8 
below.  However,  use  of  the  APEX  method  in  the  manner  to  be  described  below  is  intended 
only  for  images  where  Fourier  behavior  is  similar  to  that  shown  in  Figure  3(a). 

5.  SECB  deblurring  and  diffusion  equations 

The  SECB  method  is  a direct  (non-iterative)  FFT-based  image  deblurring  technique  designed 
for  equations  in  the  form  of  Eq.  (5),  where  h(x,y ) is  assumed  known  and  belongs  to  G. 
A complete  discussion  of  that  method,  together  with  error  bounds  and  comparisons  with 
other  methods,  may  be  found  in  Ref.  27.  Significantly,  the  SECB  method  does  not  impose 
smoothness  requirements,  such  as  prescribed  bounds  on  the  Laplacian  or  other  derivatives 
of  the  unknown  image  f(x,y).  This  is  an  important  consideration  since  many  images  have 
sharp  edges  and  other  localized  non-differentiable  features.  In  addition,  knowledge  of  the 
actual  statistical  character  of  the  data  noise  n(x , y)  in  Eq.  (6)  is  not  required,  and  the  noise 
may  be  ’multiplicative.  However,  an  estimate  of  the  L2  norm  of  n(x , y)  is  required. 

Considerable  experience  has  been  accumulated  with  the  SECB  method.  That  experience 
indicates  that  the  SECB  method  can  often  recover  fine-scale  features  in  cases  where  this 
is  not  feasible  with  iterative  methods  such  as  the  Lucy- Richardson,  Maximum  Entropy,  or 
Marquina-Osher  methods.  Documented  numerical  experiments  supporting  these  claims  may 
be  found  in  Refs.  2.  12,  and  27. 

Class  G psfs  are  the  Green's  functions  for  certain  linear  fractional  diffusion  equations.  As 
a consequence,  the  blurred  noisy  image  g{x,y)  on  the  right  of  Eq.  (5)  can  be  interpreted  as 
the  noise  corrupted  solution,  at  time  t = 1,  of  the  diffusion  initial  value  problem 

% = -IA,(-A)\  0<t<  1. 

OT  i=  1 

u{x,y,  0)  = fe{x,y ),  (13) 

where  A,  = a,-( 4ttj J — , and  A denotes  the  Laplacian.  When  the  exact  initial  value  fe{x,  y) 
is  given.  u(x . y . t)  — Hl  fe  is  the  solution  of  Eq.  (13)  at  time  t , and  u(x1  y , 1)  = ge(x , y),  in 
agreement  with  Eq.  (7).  Here.  H1  is  the  convolution  integral  operator  defined  in  Eq.  (8). 

Solving  the  deconvolution  problem  in  Eq.  (5)  is  equivalent  to  solving  the  ill-posed  back- 
wards in  time  problem  in  Eq.  (13).  namely,  given  the  noisy  data  g(x,y)  at  time  t = 1,  find 
an  approximation  f(x.y)  to  the  initial  data  fe(x,y)-  The  SECB  method  is  a regularization 
method  for  solving  that  ill-posed  diffusion  problem,  one  that  takes  into  account  the  presence 
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of  noise  in  the  blurred  image  data  g(x,  y)  at  t = 1.  The  SECB  deblurred  image  /T(x,  y)  is  an 
approximation  to  fe(x,y)  that  is  obtained  in  closed  form  in  Fourier  space.  With  z denoting 
the  complex  conjugate  of  z, 


m,v) 


h(i-v)9^-n) 

\ks,v)\2  + K~2\1  - h‘(Z,tt)\3' 


(14) 


leading  to  P(x,y)  upon  inverse  transforming.  Here,  the  positive  constants  s<  1.  and  I\. 
are  regularization  parameters,  chosen  on  the  basis  of  prior  information  as  discussed  in  Ref. 
2.  Typical  values  used  in  this  paper  might  be  s = 0.01,  K = 1000.  As  in  Eq.  (8),  we  also 
form  and  display 

u\x,y,t)  = Hf  f\x,  y),  (15) 


for  selected  decreasing  values  of  t lying  between  1 and  0.  This  simulates  marching  backwards  in 
time  in  (13),  and  allows  monitoring  the  gradual  deblurring  of  the  image.  As  t —>  0 the  partial 
restorations  u\x,y,t)  become  sharper.  Such  slow  motion  deconvolution  allows  detection  of 
features  in  the  image  before  they  become  obscured  by  noise  or  ringing  artifacts.  As  will 
be  seen  below,  such  marching  backwards  in  time  is  a vital  element  in  the  APEX  method. 
Diagnostic  statistical  information  about  uflx,  y.  t)  can  also  be  calculated  for  selected  values 
of  t as  t — » 0.  Of  particular  interest  are  the  discrete  Ll  norm,  defined  as  follows  for  2Ar  x 2N 
images 

II  0(t)  ||L1=  (2AT)-2  Y \V{x ,y,t)\,  (16) 

x,y=l 


and  the  discrete  total  variation  or  TV  norm,  which  measures  image  gradients 


2N—1 


0(t)  \\tv—  (2Ar)  2 Y ({4(x.  y.  t)}2  + y,  f)}' 

x.y=l 


A 1/2 


(17) 


where 


u'r(.r.  y.t)  = (2N)  1(u'  (x  + 1,  y,  t)  — u\x,  y,t)) 

u\(x,y,t)  = (2Ar)_1  ( u'(x,  y+  1 ,t)  -u*{x,y,t)}  (18) 

I11  blind  deconvolution  applications  of  the  SECB  method.  APEX-detected  values  for  cq,  d;. 
are  used  to  form  the  2N  x 2Ar  array  in  Eq.  (9).  This  is  input  into  Eq.  (14),  and  inverse  FFT 
algorithms  are  then  used  to  obtain  v2(x,  y.t)  in  Eq.  (15).  This  may  result  in  individual  pixel 
values  that  are  negative.  Accordingly,  all  negative  values  are  reset  to  the  value  zero.  For  such 
non-negative  image  data,  the  discrete  L 1 norm  ||  id{t)  \\w  in  Eq.  (16)  is  proportional  to  the 
total  flux.  In  a well-behaved  deconvolution  process,  this  total  flux  should  be  conserved,  and 
||  uJ(t)  \\Li  should  remain  constant,  as  t — > 0.  At  the  same  time,  the  discrete  image  gradient 
norm  ||  uflt)  \\tv  in  Eq.  (17)  should  increase  monotonically  as  t — > 0.  reflecting  the  gradual 
sharpening  of  edges  and  other  localized  singularities  in  the  restored  image. 
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6.  A priori  non- uniqueness  in  blind  deconvolution 

Blind  deconvolution  seeks  to  deblur  an  image  without  knowing  the  cause  of  the  blur.  This 
is  a difficult  mathematical  problem  in  which  severe  ill-conditioning  is  compounded  with 
non-uniqueness  of  solutions.  A-priori  constraints  can  reduce,  but  not  entirely  eliminate,  the 
multiplicity  of  solutions.  While  many  of  these  solutions  are  physically  meaningless  and  can  be 
rejected  on  physical  grounds,  there  often  remain  infinitely  many  visually  distinct,  physically 
meaningful  solutions.  Consider  the  experiment  in  Figure  4. 

The  sharp  512  x 512  Sydney  image  fe(x,  y)  in  Figure  4(a)  was  synthetically  blurred  by 
convolution  with  a Lorentzian  density  h(x,  y)  with  Qo  = 0.075,  3o  = 0.5.  This  produced 
the  blurred  image  ge{x,y)  in  Figure  4(b).  To  avoid  distractions  caused  by  noise,  the  blurred 
image  ge(x,y)  in  this  experiment  was  computed  and  stored  in  64-bit  precision.  Deblurring 
Figure  4(b)  with  the  correct  psf  values  a = 0.075,  3 - 0.5,  produces  Figure  4(c).  This  is  in 
excellent  visual  agreement  with  fe(x.  y)  in  Figure  4(a),  as  expected.  However,  Figure  4(d), 
obtained  from  Figure  4(b)  using  the  "incorrect"  psf  values  a = 0.195,/?  = 0.4,  appears  even 
sharper!  It  is  not  evident  how,  or  why , one  would  eliminate  the  reconstruction  in  4(d).  Both 
deblurred  images  were  obtained  using  the  SECB  method  with  s = 0.001  and  K = 10000. 
One  dimensional  cross  sections  of  the  two  distinct  psfs  used  in  Figure  4 are  displayed  in 
Figure  5.  These  psfs  exhibit  distinct  heavy  tail  behavior  not  shown  in  Figure  5.  The  two 
restorations  also  have  distinct  Ll  and  TV  norms,  as  shown  in  Table  1. 

Note  that  Figure  4(d)  was  obtained  using  a specific  pair  (ct,/3)  where  a > <n0,  and 
i < 30.  In  fact,  there  are  infinitely  many  other  specific  pairs  (a,/?),  capable  of  producing 
distinct,  high  quality  reconstructions  from  the  same  blurred  image  ge(x,y)  in  Figure  4(b). 
These  reconstructions  may  differ  markedly  from  one  another  at  individual  pixels,  while  being 
correct  visual  representations  of  the  object  that  was  imaged.  This  is  an  inherent,  a-priori, 
non-uniqueness  property  of  the  blind  deconvolution  problem,  independently  of  any  particular 
algorithm  that  might  be  used  to  solve  that  problem. 

This  situation  is  reminiscent  of  the  multitude  of  distinct  images  that  often  exist  for  some 
unique  astronomical  objects,  such  as  the  Whirlpool  Galaxy  (M51),  for  example.  In  that 
. these  noticeably  different  photographic  representations  of  the  identical  object  are  all 
physically  meaningful  and  visually  correct. 

I he  non- uniqueness  of  good  solutions  to  the  blind  deconvolution  problem  has  not  been 
fully  explored  in  the  literature.  When  a blind  algorithm  produces  a unique  solution,  this 
may  only  indicate  that  that  solution  is  the  only  one  accessible  to  that  particular  algorithm. 
Conceivably,  there  may  be  numerous  additional  good  solutions  that  remain  inaccessible  to 
the  algorithm.  And  some  of  these  reconstructions  may  exhibit  features  of  great  interest. 

A basic  property  of  the  APEX  method  is  that  it  generally  provides  several  psfs  that  can  be 
used  to  obtain  useful  reconstructions  of  a given  blurred  image.  As  in  the  example  above,  these 
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(c)  (d) 


Fig.  4.  Non-uniqueness  in  blind  deconvolution.  Distinct  point  spread  functions 
exist  that  produce  distinct  high  quality  reconstructions  from  the  same  blurred 
image,  (a)  Original  sharp  512  x 512  Sydney  image,  (b)  Synthetically  blurred 
Sydney  image  created  by  convolution  with  Lorentzian  density  obtained  by 
choosing  a = 0.075,  (3  = 0.5  in  Eq.  (2).  (c)  Deblurring  of  image  (b)  using 
correct  otf  parameters  a = 0.075,  (3  = 0.5.  (d)  Deblurring  of  image  (b)  using 
"incorrect"  otf  parameters  a = 0.195,  j3  = 0.4.  Deblurred  images  obtained 
using  SECB  procedure  in  Section  5,  with  s = 0.001  and  K = 10000. 
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TWO  DISTINCT  PSFS  THAT  DEBLUR  SYDNEY  IMAGE 


Fig.  5.  Two  distinct  point  spread  functions  that  deblnr  image  (b)  in  Figure  4. 
Curves  C and  D are  1-D  cross  sections  of  the  512  x 512  psfs  that  respectively 
produced  images  (c)  and  (d)  in  Figure  4.  These  psfs  also  exhibit  distinct  heavy 
tail  behavior. 


TABLE  1 

Behavior  in  deblurred  images  in  Figure  4- 


Restoration 

Q,  (3 

L1  norm, 

TV  norm 

Image  (C) 

a = 0.075.  (3  = 0.500 

173 

6419 

Image  (D) 

a = 0.195,  13  = 0.400 

171 

7500 
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reconstructions  will  differ  from  one  another  at  individual  pixels  while  being  visually  correct. 
As  is  well-known,  a-priori  knowledge  about  the  desired  solution  is  a necessary  ingredient 
for  solving  ill-posed  inverse  problems.  Such  knowledge  is  expected  to  guide  the  user  in  his 
selection  of  the  best  solution,  out  of  the  multiplicity  of  good  solutions. 

7.  Slow  motion  blind  deconvolution  and  the  APEX  method 

The  following  observations  underlie  the  APEX  method.  In  the  basic  relation 

9(x,y)  = h(x,  y)  ® fe(x,y)  + n{x,y),  (19) 

we  may  safely  assume  that  the  noise  n(x,  y)  satisfies 

/ | n(x,  y)\dxdy  < / fe(x,  y)dxdy  = 7 > 0,  (20) 

Jr2  Jr2 

so  that, 

!«*(?,  ??)|  < 1.  (21) 

Consider  the  ease  where  the  otf  is  a pure  Levy  density  /?(?,'/)  = exp{—  a(f2  + J/2)'*}-  Since 

g = ge  + n 

In  \g'(Z,v)\  = !n  I exp{-a(y  + if)s}fe  (t,v)  + »*(f,/?)|.  (22) 

Let  fi  = ( ( f . r\ ) | f2  + ip  < or}  be  a neighborhood  of  the  origin  where 

exp{—a(£,2  + i?2)'3}|/e*(f, tj)|  » |h*(Ch)l-  (23) 

Such  an  Q exists  since  Eq.  (23)  is  true  for  £ = r]  = 0 in  view  of  Eq.  (21).  The  radius  u > 0 
of  Q decreases  as  a,  /?,  and  n increase.  However,  in  many  applications,  o,  /?,  and  n(x,  y)  are 
sufficiently  small  that  Q extends  into  the  high-frequency  range.  For  (£,  //)  6 Q we  have 

~ -Vt2  + v2)i3  + 'll  I A (C-h)l-  (24) 

Because  of  the  radial  symmetry  in  the  psf.  it  is  sufficient  to  consider  Eq.  (24)  along  a single 
line  through  the  origin  in  the  (£,  rf)  plane.  Choosing  the  line  ?/  = 0.  we  have 

'n  !<?*(£.  0)|  « -a|^|23  + ln|/e*(?.0)|.  If  | < w.  (25) 

Some  type  of  a-priori  information  about  fe{x.  y)  is  necessary  for  blind  deconvolution.  In  Eq. 
(25),  knowledge  of  In  \ fe  (£,  0)|  on  |£|  < x would  immediately  yield  cx|£|2J  on  that  interval. 
Moreover,  any  other  line  through  the  origin  could  have  been  used  in  Eq.  (24).  However, 

A 

In  |/e  (£,  0)1  is  highly  oscillatory,  and  such  detailed  knowledge  is  unlikely  in  practice.  Nor 
is  it  actually  necessary.  Much  cruder  knowledge,  in  the  form  of  the  smooth  curve  T that 
best  approximates  In  \ fe  (£,0)|  in  the  least  squares  sense,  turns  out  to  be  sufficient.  Indeed, 
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knowledge  of  the  smooth  curve  T is  the  basis  for  the  BEAK  method 1 of  determining  a and 
3 from  Eq.  (25).  However,  when  T is  not  available,  the  APEX  method  must  identify  a useful 
psf  from  Eq.  (25).  using  more  elusive  information  about  In  |/e  (£,  0)|.  To  compensate  for 
this  handicap,  the  SECB  marching  backwards  in  time  option  in  Eq.  (15)  is  used,  together 
with  visual  monitoring  of  the  unfolding  deconvolution.  Accompanying  diagnostic  statistical 
information  as  t — > 0.  such  as  the  discrete  norms  ||  id(t)  ||Li,  and  ||  vf(t)  \\tv,  in  Eqs. 
(16)  and  (17),  provide  the  means  for  readjusting  initially  detected  psf  parameters,  a and  /?, 
and  enforcing  conservation  of  total  flux.  The  method  assumes  that  fe(x,y)  is  a recognizable 
object,  and  may  require  several  interactive  trials  prior  to  locating  a suitable  psf.  As  previously 
noted,  such  trial  SECB  restorations  are  easily  obtained. 

8.  Conservation  of  total  flux 

A jjc 

In  the  absence  of  the  smooth  least  squares  fit  T,  we  replace  In  \ fe  (£,  0)|  by  a negative  constant 
—A  in  Eq.  (25).  For  any  A > 0,  the  approximation 

ln|s*(e0)|  « -a|£|2<3  -A,  (26) 

is  not  valid  near  £ = 0,  since  the  curve  u(£)  = — o|£|2/3  — A,  has  —A  as  its  apex.  Choosing  a 
value  of  A > 0,  we  best  fit  In  |g*(£,0)|  with  u(£)  = —a |£|2/?  — A on  the  interval  |£|  < w, 
using  nonlinear  least  squares  algorithms.  This  is  illustrated  in  Figure  3(b).  The  resulting  fit 
is  close  only  for  f away  from  the  origin.  The  returned  values  for  a and  /3  are  then  used  in  the 
SECB  deblurring  algorithm.  Different  values  of  A return  different  pairs  (a,/3).  Experience 
indicates  that  useful  values  of  A generally  lie  in  the  interval  3 < A < 6.  Increasing 
the  value  of  A decreases  the  curvature  of  u(£)  at  £ = 0,  resulting  in  a larger  value  of 
6 together  with  a smaller  value  of  a.  A value  of  A > 0 that  returns  (3  > 1 is  clearly 
too  large,  as  8 > 1 is  impossible  for  probability  density  functions.15  Decreasing  A has  the 
opposite  effect,  producing  lower  values  of  (3  and  higher  values  of  a.  As  a rule,  deconvolution 
is  better  behaved  at  lower  values  of  /3  than  it  is  when  (3  « 1.  A significant  discovery  is  that 
an  image  blurred  with  a pair  (a0,  8q)  can  often  be  successfully  deblurred  with  an  appropriate 
pair  (a,  ,3).  where  a > a0  and  (3  < /3q  . An  example  of  this  phenomenon  was  shown  in 
Figure  4(d)  in  connection  with  the  blurred  Sydney  image.  An  effective  interactive  framework 
for  performing  the  above  least  squares  fitting  is  the  fit  command  in  DATAPLOT28  This  is  a 
high-level  English-syntax  graphics  and  analysis  software  package  developed  at  the  National 
Institute  of  Standards  and  Technology.  This  software  tool  was  used  throughout  this  paper. 

The  following  version  of  the  APEX  method  has  been  found  useful  in  a variety  of  image 
enhancement  problems  where  the  image  g(x,  y)  is  such  that  In  | <?*(£,  0)|  is  generally  globally 
monotone  decreasing  and  convex,  as  shown  in  Figure  3(a).  Choose  a value  of  A > 3 in  Eq. 
(26),  so  that  the  least  squares  fit  develops  a well-formed  cusp  at  £ = 0,  as  shown  in  Figure 
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3(b).  Using  the  returned  pair  (a,  (3)  in  the  SECB  method,  obtain  a sequence  u^(x.yA)  of 
partial  restorations  as  in  Eq.  (15),  as  t decreases  from  t — 1.  With  a good  choice  of  A , the 
total  flux  norm,  ||  uJ(t)  |fti,  should  remain  constant  or  increase  very  slowly  as  t decreases, 
while  the  image  gradient  norm,  ||  vf{t)  \\tv,  should  increase  monotonically  as  t decreases 
from  t — 1. 

Most  often,  the  initially  detected  value  of  a turns  out  to  be  too  large.  The  corresponding  psf 
is  then  too  wide  in  physical  (x,  y)  space,  or,  equivalently,  the  otf  is  too  narrow  in  Fourier  (£,  y) 
space.  Theoretically,  use  of  too  wide  a psf  all  the  way  to  t = 0,  implies  sharpening  features 
that  may  have  already  become  infinitely  sharp  at  some  tG  > 0.  In  practice,  this  leads  to  severe 
ringing  and  other  undesirable  artifacts  at  t = 0.  indicating  that  continuation  backwards  in 
time  has  proceeded  too  far.  An  accompanying  symptom  of  this  ill-behaved  deconvolution,  is 
that  the  total  flux  norm  ||  uJ(t)  |fti  does  not  remain  constant,  but  increases  appreciably  as 
t — » 0.  Choosing  a new  and  larger  value  of  A in  Eq.  (26),  returns  a smaller  a,  but  with  a 
larger  [3.  A useful  strategy  is  to  locate  a pair  (<n,/5)  such  that  ||  v3(t)  ||Li  increases  slowly 
enough  as  t decreases,  that  its  value  at  t = ta  = 0.5,  say,  is  only  a very  few  percent  more 
than  its  initial  value  at  t — 1.  In  that  case,  the  deconvolution  is  terminated  at  t — ta.  To 
enforce  total  flux  conservation,  the  resulting  image  at  ta  is  rescaled  by  multiplying  it  by  the 
constant  Ca  =||  iT(  1)  ||^i  / ||  iT(ta)  \\Li.  Ideally,  Ca  should  be  very  close  to  unity.* 

Marching  backwards  in  time  allows  for  simultaneous  sampling  of  numerous  values  of  a 
while  keeping  f3  fixed.  Terminating  the  deconvolution  at  t — ta  > 0,  is  equivalent  to  read- 
justing the  original  a while  keeping  the  same  value  of  (3.  If  the  pair  (a,  j3)  produces  a high 
quality  restoration  at  t = ta  > 0,  the  pair  (q*,/5),  where  a*  = (1  — ta) a,  will  produce  the 
same  quality  results  at  t — 0.  We  therefore  distinguish  between  the  originally  detected  a, 
and  the  effective  a , a *.  In  general,  there  will  be  many  values  of  A in  (26)  returning  pairs 
(ct,/3)  that  produce  good  reconstructions  at  some  tay  > 0.  A large  number  of  distinct  pairs 
(<n*,/3)  can  thus  be  found  that  produce  useful,  but  distinct,  results  at  t = 0. 

Ideally,  successful  APEX  blind  deconvolution  should  incorporate  three  elements:  clear 
visual  evidence  of  sharpening,  accompanied  by  a substantial  increase  in  TV  norm,  and 
conservation  of  L 1 norm. 

We  have  been  assuming  /?(£,  ??)  to  be  a pure  Levy  otf  in  Eq.  (19).  The  procedure  is 
very  similar  for  the  more  general  class  G otfs  in  Eq.  (4).  Here,  given  prior  starting  values 
for  the  cq,  ft,  i = 1,  J,  we  best-fit  ln|p*(£,0)|  with  — XjLi  cp|£|2/?z  — A,  , with  suitably 
preselected  A > 3.  This  returns  J initially  detected  pairs  (cq,ft).  As  before,  by  monitoring 
the  deconvolution  process  and  terminating  it  at  the  appropriate  time  ta  > 0.  we  arrive  at 
effective  values  a*  = (1  — ta)al,  such  that  the  J pairs  (a-, ft)  produce  useful  sharpening  at 

'It  is  occasionally  beneficial  to  allow  more  aggressive  deblurring,  with  the  L 1 norm  increasing  by  as  much 
as  10%  prior  to  rescaling,  in  order  to  bring  out  important  fine  scale  structural  details. 
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t = 0.  It  should  be  noted  that  in  most  applications  of  the  APEX  method  considered  to  date, 
including  those  in  the  present  paper,  high  quality  reconstructions  were  obtained  using  the 
simplest  version  of  that  method,  where  J = 1.  This  indicates  that  in  many  applications,  a 
single  pure  Levy  stable  otf  can  often  be  found  that  sufficiently  well- approximates  the  system’s 
more  complex  composite  otf. 

All  point  spread  and  optical  transfer  functions  depicted  in  this  paper,  including  those  in 
Figure  5.  are  based  on  effective  Levy  parameter  values  (a*,/?),  producing  optimal  recon- 
structions at  t = 0. 

9.  Applications  to  gray  scale  galaxy  images 

Our  first  example,  in  Figure  6(a),  is  a 1024  x 1024  8-bit  gray  scale  image  g(x,y)  of 
the  spiral  galaxy  M101.  This  is  adapted  from  a similar  size  color  JPEG  image  obtained 
by  Jacoby.  Bohannan,  and  Hanna,  Kitt  Peak  National  Optical  Astronomy  Observatory, 
(NOAO/AURA/NSF).  A plot  of  ln|g*(£,0)|  was  shown  earlier  in  Figure  3(a).  Using 
A = 3.85.  we  best-fit  ln|g*(£,0)|  with  — a|£|2i  — A on  |f|  < 500.  The  fit  develops  a 
well-formed  cusp  at  £ = 0,  as  shown  in  Figure  3(b),  and  returns  a = 0.385,  [3  = 0.206. 
The  patterns  illustrated  in  Figure  3 are  typical  of  all  the  images  discussed  in  this  paper. 
Using  these  Levy  parameters  in  the  SECB  method  with  s = 0.01,  K = 1300,  we  find  that 
at  first,  ||  i A(t)  || £i  increases  slowly  as  t decreases,  from  an  initial  value  of  12.84  at  t = 1, 
to  a value  of  12.96  at  t — 0.65.  Thereafter,  ||  vf(t)  ||Li  increases  more  rapidly.  At  the  same 
time,  ||  uT(t)  || tv  increases  monotonically  from  2134  at  t = 1,  to  6747  at  t = 0.65,  i.e.,  a 
threefold  increase  in  gradient  norm.  Deconvolution  was  terminated  at  ta  = 0.65,  and  the 
effective  value  of  a is  a*  = 0.135.  The  APEX-processed  image,  shown  in  Figure  6(b),  was 
rescaled  so  as  to  have  the  same  L 1 norm  as  Figure  6(a). 

Our  second  example,  in  Figure  7(a),  is  a 1024x1024  8-bit  gray  scale  image  of  the  spiral 
galaxy  M51.  This  is  adapted  from  a similar  size  color  JPEG  image  obtained  by  Rector  and 
Ramirez.  Kitt  Peak  National  Optical  Astronomy  Observatory,  (NOAO/AURA/NSF).  Here, 
there  is  very  substantial  documented  APEX  sharpening,  and  the  deconvolved  image  in  Figure 
7(b)  very  visibly  improves  on  the  original.  With  A=5.0,  least  squares  fitting  on  |£|  < 500, 
returned  a = 0.364.  8 = 0.218.  This  was  input  into  the  SECB  method  with  s=0.01,  K=1300. 
Deconvolution  was  terminated  at  ta  = 0.48,  leading  to  an  effective  a*  = 0.189.  Total  flux 
||  vJ(t)  \\li  increased  very  slightly,  from  30.58  at  t — 1,  to  30.82  at  ta  = 0.48.  However, 
there  was  a corresponding  eightfold  increase  in  ||  iA(t)  \\tv,  from  1948  at  t = 1,  to  16516  at 
ta  = 0.48.  Both  images  in  Figure  7 have  identical  V norms. 

Our  next  example,  in  Figure  8(a).  is  a 1024  x 1024  8-bit  gray  scale  image  of  the  spiral 
galaxy  M74.  This  is  adapted  from  a JPEG  color  image  taken  in  August  2001  by  the  GMOS 
Team  at  the  Gemini  Observatory.  Manna  Kea,  Hawaii.  With  A = 4.25,  least  squares  fitting 
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Fig.  6.  APEX  blind  deconvolution  of  M101  image,  (a)  Original  1024  x 1024 
M101  image,  obtained  by  Jacoby,  Bohannan,  and  Hanna,  Kitt  Peak  National 
Observatory,  (NOAO/AURA/NSF).  (b)  APEX  processed  image  is  noticeably 
sharper.  Both  images  have  identical  AS  ‘total  flux'  norms,  but  TV  'gradient' 
norm  in  image  (b)  is  three  times  larger  than  in  image  (a). 


Fig.  7.  APEX  blind  deconvolution  of  M51  image,  (a)  Original  1024  x 1024 
M51  image  obtained  by  Rector  and  Ramirez,  Kitt  Peak  National  Observatory, 
(NOAO/AURA/NSF).  fb)  APEX  processing  very  significantly  improves  orig- 
inal. 'Total  flux'  L]  norms  of  images  (a)  and  (b)  are  equal,  but  ‘gradient/  TV 
norm  in  (b)  is  more  than  eight  times  nftger  than  in  (a). 


Fig.  8.  APEX  blind  deconvolution  of  M74  image,  (a)  Original  1024  x 1024  M74 
image  obtained  by  GMOS  Team  at  Gemini  Observatory,  Mauna  Ivea.  Hawaii, 
(b)  APEX  processing  significantly  improves  original.  ‘Total  flux'  L1  norms  in 
images  (a)  and  (b)  are  equal,  but  ‘gradient-  TV  norm  in  (b)  is  four  times  larger 
than  in  (a). 
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Fig.  9.  APEX  processing  leads  to  significant  change  in  high  frequency  Fourier 
behavior,  (a)  Before  and  after  for  M51  image,  (b)  Before  and  after  for  M74 
image. 


of  In  |<7*(£.0)|  with  — a|£|Jl  — A,  on  |£|  < 500,  returned  a = 0.857,  ft  = 0.157.  Here,  more 
aggressive  deblurring  was  permitted  prior  to  termination.  With  s = 0.01  and  K = 500  in 
the  SECB  method,  ||  id(t)  ||^i  increased  by  7%,  from  61.22  to  65.68,  prior  to  termination 
at  ta  = 0.65.  The  effective  value  of  a is  a*  = 0.3,  and  there  was  a corresponding  fourfold 
increase  in  ||  uJ(t)  \\tv,  from  4670  to  20173.  The  APEX-processed  image  in  Figure  8(b)  was 
rescaled  so  as  to  have  the  same  Ll  norm  as  Figure  8(a). 

In  Figure  9,  In  | <?*(£,  0)|  is  plotted  on  |£|  < 500,  for  the  M51  and  M74  images,  before  and 
after  APEX  processing.  Evidently,  APEX  processing  amplifies  high  frequency  components 
quite  significantly.  This  amplification  is  carefully  orchestrated,  takes  place  in  a stable,  co- 
herent fashion,  and  enables  recovery  of  the  delicate  fine  structures  and  other  features  that 
are  evident  in  Figures  7(b)  and  8(b).  These  before  and  after  Fourier  patterns  are  typical  of 
all  the  images  shown  in  this  paper. 

10.  APEX  processing  of  color  imagery 

Blind  deconvolution  of  color  imagery  is  a subject  that  is  still  very  much  in  its  infancy.  Major 
difficulties  arise  from  the  need  to  identify  the  distinct  point  spread  functions  associated  with 
each  color  component.  More  serious  difficulties  arise  from  the  possibility  of  unbalanced  blind 
sharpening  of  individual  color  components.  Conceivably,  after  a long  and  uncertain  iterative 
process,  the  reconstituted  color  image  may  turn  out  to  exhibit  physically  false  colors,  such 
as  a green  sky.  or  a purple  sea.  A fruitful  mathematical  framework  wherein  the  blind  color 
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problem  can  be  effectively  tackled,  has  not  yet  been  formulated. 

One  approach  to  color  image  processing  traces  its  origin  to  high  energy  physics  and 
string  theory.29-31  ffere,  a color  image  is  viewed  as  a 2D  manifold  in  5 D space,  namely, 
{x,y.  R(x,  y),G(x,  y).  B(x.  y)}.  where  i?,  G.  B are  the  red,  green,  and  blue  components 
of  the  color  image  g(x,y).  The  so-called  Polyakov  functional  is  then  defined  on  this  man- 
ifold, and  gradient  descent  minimization  of  this  functional  is  implemented.  This  leads  to 
the  Beltrami  flow  equations,  a coupled  system  of  evolutionary  nonlinear  partial  differential 
equations  for  the  three  time- dependent  images  R(x,  y.t),  G(x,y.t ),  B(x,  y,  t).  That  system 
is  then  solved  forward  in  time  numerically,  until  a steady-state  is  reached.  This  formalism 
has  been  applied  successfully  to  color  image  denoising.  With  considerable  skill,  such  an  ap- 
proach might  possibly  be  elaborated  into  a blind  deconvolution  procedure.  However,  the 
computational  effort  required  to  process  large  size  imagery  would  be  challenging. 

A remarkable  property  of  the  APEX  method  is  the  ease  with  which  it  can  be  applied  to 
color  imagery,  and  the  plausibility  of  the  ensuing  results.  Clearly,  the  ability  to  try  numerous 
parameter  values  in  quasi  real-time  is  of  vital  significance.  Indeed,  efficient  exploration  in 
parameter  space  is  often  the  key  to  the  successful  solution  of  ill-posed  inverse  problems. 

The  most  natural  way  to  use  the  APEX  method  is  to  first  decompose  the  blurred  color 
image  into  its  three  RGB  components,  apply  the  method  to  each  component  in  turn,  and  then 
reconstitute  the  deblurred  image.  For  each  RGB  component  , visual  monitoring  of  the  partial 
deconvolution  uflx,  yfl)  in  Eq.  (15)  as  t —>  0,  is  accompanied  by  the  calculated  diagnostic 
quantities  ||  uflt)  ||Li  and  ||  uflt)  \\tv-  As  in  the  case  of  gray  scale  imagery  discussed  above, 
total  flux  conservation  in  each  RGB  component  is  enforced  by  terminating  deconvolution  at 
some  appropriate  time  ta  > 0,  and  rescaling  the  image  by  multiplication  by  the  constant 
Ca  =||  td(  1)  ||Li  / ||  iflfla)  || l1-  In  this  way,  individual  Levy  pairs  (a*,/?)  are  detected  for 
each  RGB  component,  often  leading  to  distinct  otfs  for  each  color.  This  methodology  has 
also  been  found  to  maintain  the  balance  of  colors  in  all  of  the  many  examples  to  which  it 
has  been  applied.  We  shall  now  demonstrate  this  on  several  color  images,  including  some 
spectacular  Hubble  Space  Telescope  images. 

Our  first  color  image,  in  Figure  10(a),  is  an  Ektachrome  1024  x 727  tiff  image  of  the 
Southern  Pinwheel  galaxy  M83.  This  was  downsized  from  an  original  3500  x 2485  tiff 
image  taken  by  Bill  Schoening  at  Ivitt  Peak  National  Optical  Astronomy  Observatory, 
(NOAO/AURA/NSF).  After  decomposition  into  RGB  components,  APEX  least  squares  fit- 
ting on  |£|  < 500,  with  A = 3.85,  was  applied  to  each  component.  The  returned  values  for 
a and  3 were  then  input  into  the  SECB  method  with  s — 0.01.  K = 1300.  For  the  red 
component,  a — 0.356.  3 = 0.195,  and  ||  uJ(t)  ||Li  increased  by  about  2%,  from  15.02  to 
15.36,  prior  to  termination  at  ta  = 0.65.  The  effective  a in  this  case  is  a*  = 0.124.  There 
was  a corresponding  threefold  increase  in  ||  uflt)  \\tv-,  from  2541  to  7698.  For  the  green 
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Fig.  10.  APEX  processing  significantly  sharpens  Southern  Pinwheel  M83  Ek- 
taehrome  image.  Original  (a)  was  obtained  by  Bill  Schoening,  Ivitt  Peak  Na- 
tional Observatory,  (XOAO/AURA/NSF).  Both  images  have  equal  L 1 ‘total 
flux'  norms  in  each  RGB  component,  but  component  TV  ‘gradient'  norms  in 

enhanced  image  (b)  are  three  times  larger  than  in  (a). 
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Fig.  11.  APEX  processing  enhances  Hubble  Space  Telescope  image  of 
NGC2207  involving  two  merging  galaxies.  Original  (a)  was  obtained  by  NASA, 
ESA,  and  the  Hubble  Heritage  Team  (STSci/AURA).  Both  images  have  equal 
L 1 ‘total  flux"  norms  in  each  RGB  component,  but  component  TV  ‘gradient' 
norms  in  enhanced  image  (b)  are  almost  three  times  larger  than  in  (a). 
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Fig.  12.  APEX  blind  deconvolution  enhances  Hubble  Space  Telescope  image  of 

Orion  Reflection  Nebula.  XGC  1999.  Original  (a)  was  obtained  by  NASA  and 

the  Hubble  Heritage  Team  (STSci/AURA).  Both  images  have  equal  V ‘total 
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flux  norms  in  each  RGB  component,  but  component  TV  ‘gradient’  norms  in 
enhanced  image  (b)  are  5.6  times  larger  than  in  (a). 


component,  a = 0.506,  ,3  = 0.171,  and  the  Ll  norm  increased  by  5%,  from  23.76  to  24.99. 
prior  to  termination  at  tG  = 0.65.  Here,  a*  — 0.177.  There  was  again  a threefold  increase  in 
TV  norm,  from  4393  to  14613.  For  the  blue  image,  a = 0.571,  f3  — 0.161.  and  ||  uJ(t)  ||Li 
increased  by  6.5%,  from  30.60  to  32.59,  prior  to  termination  at  tG  = 0.65.  This  gives  a*  = 0.2. 
Once  again,  there  was  a threefold  increase  in  ||  uJ(t)  || tv-,  from  3606  to  11778. 

In  this  example,  the  green  image  otf  almost  coincides  with  the  blue  image  otf,  and  both  lie 
below  the  red  image  otf.  Thus,  APEX  methodology  perceived  the  red  component  to  be  less 
blurred  than  the  other  two  components,  and  it  processed  the  image  accordingly.  All  three 
components  were  rescaled  to  preserve  V norms,  prior  to  reconstitution  into  Figure  10(b). 

The  next  example  is  a Hubble  Space  Telescope  image  of  NGC2207,  involving  two  merging 
galaxies.  That  image  forms  part  of  the  Hubble  Heritage  Gallery.  The  original  full  resolu- 
tion 2907  x 1486  tiff  image  was  obtained  by  NASA,  ESA,  and  the  Hubble  Heritage  Team 
(STSci/AURA),  using  the  Wide  Field  and  Planetary  Camera  2,  (WFPC2).  This  was  stepped 
down  to  the  1024  x 523  shown  in  Figure  11(a).  We  used  A = 4.75  with  5 = 0.01,  K — 1300 
in  the  SECB  method,  and  terminated  the  process  at  tG  = 0.65  in  each  of  the  three  com- 
ponents. Here,  APEX  perceived  the  red  component  to  be  more  blurred  than  the  other  two 
components.  For  the  red  image,  a*  = 0.111,  (3  = 0.203,  and  ||  rd(£)  ||Li  increased  by  10.25%, 
from  19.11  to  21.06.  while  ||  id( t ) \\tv  increased  from  3862  to  11182,  a factor  of  2.9.  For  the 
green  image,  a * = 0.088,  f3  — 0.217,  ||  id(£)  ||Li  increased  from  17.71  to  18.82,  (6.8%),  while 
II  || TV  increased  by  a factor  of  2.8,  from  3937  to  11019.  The  blue  image  was  perceived 

to  be  the  least  blurred.  Here,  q*  = 0.052,  (3  = 0.247,  ||  id(£)  ||Li  increased  by  8.9%,  from 
14.67  to  15.97.  while  ||  vf(t)  \\tv  increased  by  a factor  of  2.4,  from  7783  to  17726.  All  three 
RGB  components  were  rescaled  to  preserve  L 1 ‘total  flux1  norms,  prior  to  reconstitution  into 
Figure  11(b). 

Our  third  example  is  again  a Hubble  Heritage  Gallery  image,  featuring  the  reflection  nebula 
in  Orion,  NGC  1999.  The  original  750  x 750  tiff  image  was  obtained  by  NASA  and  the  Hubble 
Heritage  Team  (STSci/AURA),  using  the  WFPC2  camera.  Here,  this  was  stepped  down  to 
the  512  x 512  image  shown  in  Figure  12(a).  With  A = 5.5  and  5 = 0.01,  K = 1300  in 
SECB,  deconvolution  was  unusually  well-behaved  and  uniform For  each  RGB  component. 
II  M(£)  || l1  was  very  nearly  conserved  prior  to  termination  at  t G — 0.6.  This  norm  was  near 
97  for  the  blue  image,  and  near  67  for  the  red  and  green  images.  Moreover,  all  three  otfs 
coincided  in  this  case,  as  detected  Levy  pairs  for  each  component  were  all  very  nearly  equal 
to  a*  = 0.3,  (3  = 0.17.  Again,  for  each  RGB  component,  ||  iT{i)  \\tv  increased  by  the  same 
factor  of  5.6.  from  1498  to  8539  for  red,  from  1450  to  8166  for  green,  and  from  2291  to  12890 
for  blue.  The  enhanced  image  is  shown  in  Figure  12(b). 

As  was  the  case  with  gray  scale  galaxy  images,  the  striking  improvements  in  visual  quality 
in  Figures  10(b),  11(b).  and  12(b),  appear  to  correlate  well  with  substantial  increases  in  TV 
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norms. 


11.  Advanced  Camera  for  Surveys  (ACS)  imagery 

The  WFPC2  is  Hubble's  main  camera  and  workhorse  instrument.  Our  final  two  examples 
feature  images  taken  with  the  Advanced  Camera  for  Surveys  (ACS).  That  instrument  out- 
performs all  previous  cameras  aboard  the  Hubble  Space  Telescope.  To  celebrate  Hubble’s 
fifteenth  birthday  on  April  25  2005,  NASA  released  the  sharpest-ever  color  image  of  the 
Whirlpool  Galaxy  M51.  That  image  was  recorded  with  the  ACS  Camera  by  NASA,  ESA, 
S.  Beckwith  (STScI),  and  the  Hubble  Heritage  Team  (STScI/AURA).  The  original  full  res- 
olution 7965  x 11477  tiff  image  was  stepped  down  to  the  710  x 1024  tiff  image  shown  in 
Figure  13(a).  After  decomposing  that  image  into  RGB  components,  APEX  processing  using 
A — 5.25  was  applied  to  each  component  in  turn,  with  .s  = 0.01,  K = 1300  in  the  SECB 
method.  All  three  components  behaved  very  similarly,  and  deconvolution  was  terminated  at 
ta  = 0.65  in  all  three  cases.  For  the  red  component,  a*  = 0.175,  ft  = 0.173,  and  ||  ||/,i 

increased  by  5.3%  prior  to  termination,  from  42.46  to  42.72.  However,  ||  u\t)  \\tv  increased 
by  a factor  of  3.7,  from  5170  to  19247.  For  the  green  component  a*  = 0.177,  (3  = 0.171,  and 
the  Ll  norm  increased  from  41.63  to  44.05,  a 5.8%  increase.  The  TV  norm  increased  from 
4361  to  17801.  a fourfold  increase.  For  the  blue  component,  a*  = 0.160,  f3  = 0.186,  and  the 
Ll  norm  increased  from  41.11  to  43.28,  a 5.3%  increase.  There  was  again  a fourfold  increase 
in  the  TV  norm,  from  4805  to  19839.  All  three  component  otfs  coincided  in  this  case.  Indi- 
vidual RGB  components  were  rescaled  so  as  to  preserve  L 1 norms,  prior  to  reconstitution  as 
the  APEX  image  shown  in  Figure  13(b). 

Our  last  example  involves  an  image  of  the  Tadpole  galaxy  UGC10214,  set  against  a back- 
drop that  is  said  to  contain  a "Whitman's  Sampler  of  galaxies  stretching  back  to  the  begin- 
ning of  time" . The  full  resolution  3806  x 4160  tiff  image  was  taken  with  the  ACS  Camera  by 
NASA.  STScI.  ESA,  and  the  ACS  Science  Team.  This  was  stepped  down  to  the  937  x 1024  tiff 
image  shown  in  Figure  14(a).  After  decomposing  that  image  into  RGB  components,  APEX 
processing  using  A = 5.25  was  applied  to  each  component  in  turn,  with  s = 0.01,  K = 1300 
in  the  SECB  method.  Deconvolution  was  uniformly  well-behaved,  and  was  terminated  at 
ta  = 0.675  in  all  three  components.  For  the  red  component,  q*  = 0.066,  (3  = 0.242,  and 
| n'(t)  ||^i  increased  by  2.2%  prior  to  termination,  from  23.54  to  24.05.  At  the  same  time, 
there  was  a threefold  increase  in  ||  id(£)  || TV,  from  6085  to  18606.  For  the  green  component 
a*  = 0.068.  [3  = 0.234,  and  ||  v3(t)  ||Li  increased  by  3%,  from  23.025  to  23.72.  Again,  there 
was  a near  threefold  increase  in  ||  v3(t)  \\tv , from  6970  to  19645.  For  the  blue  component, 
q*  — 0.103.  (3  = 0.201,  the  V norm  increased  by  3.2%,  from  25.25  to  26.07,  while  the  TV 
norm  increased  threefold,  from  7731  to  22756.  Again,  all  three  component  otfs  coincided.  In- 
dividual RGB  components  were  rescaled  so  as  to  preserve  L1  norms,  prior  to  reconstitution 
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Fig.  13.  APEX  processing  significantly  sharpens  15th  anniversary  Hubble 
Space  Telescope  Whirlpool  galaxy  image,  released  on  April  25  2005.  Origi- 
nal (a)  recorded  with  AOS  Camera  by  NASA,  ESA,  S.  Beckwith  (STScI), 
and  Hubble  Heritage  Team  (STScI/AURA).  Both  images  have  equal  Ll  'total 
flux'  norms  in  each  RGB  component,  but  component  TV  ‘gradient’  norms  in 
enhanced  image  (b)  are  four  times  larger  than  in  (a). 
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(b) 


Fig.  14.  APEX  processing  enhances  Hubble  Space  Telescope  Tadpole  galaxy 
image.  Original  ACS  image  (a)  taken  by  NASA,  STScI,  ESA,  and  the  ACS 
Science  Team.  Both  images  have  equal  L]  ‘total  flux'  norms  in  each  RGB 
component,  but  component  TV  ‘gracflflnt’  norms  in  enhanced  image  (b)  are 
three  times  larger  than  in  (a). 


(a) 


Fig.  15.  Extent  of  sharpening  in  APEX  processed  image  becomes  more  evident 
when  zooming  on  selected  parts  of  images  in  Figure  14.  Foreground  objects  as 
well  as  background  galaxies  in  original  (a),  are  brought  into  sharper  focus  in 
APEX  image  (b). 
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DISCRETE  FREQUENCY  (INTEGER) 


DISCRETE  FREQUENCY  (INTEGER) 


(a) 


(b) 


Fig.  16.  ID  cross  sections  of  optical  transfer  functions  that  deblurred  images 
discussed  in  sections  9 through  11.  (a)  Non-Hubble  otfs.  (b)  Hubble  otfs. 


TABLE  2 

Summary  of  APEX  experiments  in  sections  9 through  11. 


Image 

Size 

A 

ta 

a* 

0 

x TV 

M101 

1024  x 1024 

3.85 

0.65 

0.135 

0.206 

x 3 

M51  (KP) 

1024  x 1024 

5.00 

0.48 

0.189 

0.218 

x 8 

M74 

1024  x 1024 

4.25 

0.65 

0.300 

0.157 

x 4 

M83 

1024  x 727 

3.85 

0.65 

0.200 

0.161 

x 3 

Merging 

1024  x 523 

4.75 

0.65 

0.111 

0.203 

x 3 

Orion 

512  x 512 

5.50 

0.60 

0.300 

0.168 

x 6 

M51  (HST) 

710  x 1024 

5.25 

0.65 

0.177 

0.171 

x 4 

Tadpole 

937  x 1024 

5.25 

0.65 

0.068 

0.234 

x 3 
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as  the  APEX  image  shown  in  Figure  14(b). 

Zooming  on  selected  parts  of  the  images  in  Figure  14  provides  a useful  comparison,  as 
shown  in  Figure  15.  The  extent  of  sharpening  in  the  APEX  processed  image  15(b)  becomes 
clearly  evident  as  foreground  objects,  as  well  as  background  galaxies,  are  brought  into  sharper 
focus. 

The  ability  of  the  APEX  method  to  enhance  ACS  images  is  remarkable  and  unanticipated. 
Figure  16  shows  ID  cross-sectional  plots  of  the  optical  transfer  functions  that  were  detected 
and  used  to  process  all  of  the  images  discussed  in  sections  9 through  11.  Xon-Hubble  otfs 
are  shown  in  Figure  16(a)  and  Hubble  otfs  in  Figure  16(b).  The  otfs  shown  for  color  images 
are  the  ones  associated  with  their  most  blurred  RGB  component.  These  otfs  plots  are  based 
on  effective  values  (q*,/3),  that  produce  high  quality  SECB  reconstructions  at  t = 0.  These 
are  the  values  shown  in  Table  2,  which  summarizes  the  results  of  the  APEX  experiments 
described  in  sections  9 through  11.  The  last  column  in  Table  2 indicates  the  resulting 
multiplying  factors  for  TV  norm  increases.  With  an  average  value  of  i3  less  than  0.2,  the 
above  8 otfs  are  very  far  from  Lorentzian  (/ 3 = 0.5),  let  alone  Gaussian  (/?  = 1.0). 
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